In [1]:
%matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse as sp
from IPython.html import widgets # Widget definitions
from IPython.display import display, clear_output, HTML # Used to display widgets in the notebook
from IPython.html.widgets import interact, interactive, fixed
In [2]:
from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import TRES, HDF5Interface
myDataSpectrum = DataSpectrum.open("../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([22]))
myInstrument = TRES()
#myHDF5Interface = HDF5Interface("../libraries/PHOENIX_objgrid6000.hdf5")
myHDF5Interface = HDF5Interface("../libraries/PHOENIX_objgrid4000.hdf5")
myModel = Model(myDataSpectrum, myInstrument, myHDF5Interface, stellar_tuple=("temp", "logg", "Z", "vsini", "vz", "logOmega"),
cheb_tuple=("c1", "c2", "c3"), cov_tuple=("sigAmp", "logAmp", "l"), region_tuple=("h", "loga", "mu", "sigma"))
myOrderModel = myModel.OrderModels[0]
In [5]:
from StellarSpectra.model import ModelHighAccuracy
#myHDF5InterfaceHighAccuracy = HDF5Interface("../libraries/PHOENIX_submaster.hdf5")
myHDF5InterfaceHighAccuracy = HDF5Interface("../libraries/PHOENIX_LkCa15.hdf5")
#myHDF5Interface = HDF5Interface("../libraries/PHOENIX_M.hdf5")
myModelHA = ModelHighAccuracy(myDataSpectrum, myInstrument, myHDF5InterfaceHighAccuracy, stellar_tuple=("temp", "logg", "Z", "vsini", "vz", "logOmega"),
cheb_tuple=("c1", "c2", "c3"), cov_tuple=("sigAmp", "logAmp", "l"), region_tuple=("h", "loga", "mu", "sigma"))
myOrderModelHA = myModelHA.OrderModels[0]
In [3]:
myHDF5Interface.bounds
Out[3]:
In [6]:
spec = myModel.get_data()
wl = spec.wls[0]
fl = spec.fls[0]
In [7]:
def plot_diff_spec(wl, spec0, spec1, spec2=None, norm=True, abs_err=False):
'''
Plot the difference between two spectra, if normalize, divide each by their mean.
'''
if norm:
spec0 = spec0/np.mean(spec0)
spec1 = spec1/np.mean(spec1)
if np.all(spec2) is not None:
spec2 = spec2/np.mean(spec2)
fig, ax = plt.subplots(nrows=2, figsize=(11,8), sharex=True)
ax[0].plot(wl, spec0, "b", label="spec0")
ax[0].plot(wl, spec1, "r", label="spec1")
if np.any(spec2):
ax[0].plot(wl, spec2, "g", label="spec2")
ax[0].legend()
ax[0].set_ylabel("spectra")
resid1 = spec0 - spec1
if abs_err:
resid1 = np.abs(resid1/spec0)
ax[1].plot(wl, resid1 * 100, "k", label="spec0 - spec1")
if np.any(spec2):
resid2 = spec0 - spec2
if abs_err:
resid2 = np.abs(resid2/spec0)
ax[1].plot(wl, resid2 * 100, "b", label="spec0 - spec2")
ax[1].legend()
ax[1].set_xlabel(r"$\lambda$\AA")
ax[1].set_ylabel("residuals (\%)")
plt.show()
In [8]:
def create_resid_spec(spec0, spec1):
'''
Given a base spectrum and a spectrum slightly offset, return the residual spectrum between the two
'''
#Normalize both spectra
spec0 = spec0/np.mean(spec0)
spec1 = spec1/np.mean(spec1)
resid = np.abs((spec0 - spec1)/spec0)
return resid
In [9]:
def create_min_spec(base, temp_spec, logg_spec, Z_spec):
'''
Given a base spectrum and three perturbing spectra, determine for each pixel the minimum offset so that we can use
this as an error envelope in the approximate spectra.
'''
#Normalize all three spectra
base = base/np.mean(base)
temp_spec = temp_spec/np.mean(temp_spec)
logg_spec = logg_spec/np.mean(logg_spec)
Z_spec = temp_spec/np.mean(Z_spec)
#Calculate the residual spectrum for each
resid_temp = np.abs((base - temp_spec)/base)
resid_logg = np.abs((base - logg_spec)/base)
resid_Z = np.abs((base - Z_spec)/base)
#For each pixel, take the smallest value
#Vstack the arrays and take the min along axis=1
arr = np.vstack((resid_temp, resid_logg, resid_Z))
resid_min = np.min(arr, axis=0)
return resid_min
In [10]:
myModelHA.update_Model({"temp":4000, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.})
base = myOrderModelHA.get_spectrum()
myModelHA.update_Model({"temp":4010, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.})
temp_spec = myOrderModelHA.get_spectrum()
myModelHA.update_Model({"temp":4000, "logg":4.05, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.})
logg_spec = myOrderModelHA.get_spectrum()
myModelHA.update_Model({"temp":4000, "logg":4.00, "Z":-0.05, "vsini":0., "vz":0, "logOmega":0.})
Z_spec = myOrderModelHA.get_spectrum()
In [12]:
min_spec = create_min_spec(base, temp_spec, logg_spec, Z_spec)
In [13]:
plt.plot(wl, min_spec * 100)
plt.show()
In [14]:
params = {"temp":4000, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux0 = myOrderModel.get_spectrum()
In [15]:
acc_spec = create_resid_spec(base, model_flux0)
In [16]:
fig, ax = plt.subplots(nrows=2, figsize=(11,8), sharex=True)
ax[0].plot(wl, base, "b", label="base")
ax[0].plot(wl, model_flux0, "r", label="approx")
ax[0].legend()
ax[1].plot(wl, acc_spec, "r", label="approx")
ax[1].plot(wl, min_spec, "k", label="errenvelope")
ax[1].legend()
plt.show()
In [7]:
params = {"temp":6010, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux1 = myOrderModel.get_spectrum()
In [8]:
params = {"temp":6000, "logg":4.0, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModelHA.update_Model(params)
model_flux0HA = myOrderModelHA.get_spectrum()
In [31]:
plot_diff_spec(wl, model_flux0, model_flux0HA, abs_err=True)
In [24]:
plot_diff_spec(wl, base, model_flux0, abs_err=True)
In [ ]:
In [124]:
params = {"temp":2300, "logg":4.05, "Z":0.0, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux2 = myOrderModel.get_spectrum()
In [127]:
plot_diff_spec(wl, model_flux0, model_flux2, abs_err=True)
In [126]:
plot_diff_spec(wl, model_flux0, model_flux1, model_flux2, abs_err=True)
In [128]:
params = {"temp":2300, "logg":4.00, "Z":-0.05, "vsini":0., "vz":0, "logOmega":0.}
myModel.update_Model(params)
model_flux3 = myOrderModel.get_spectrum()
In [130]:
plot_diff_spec(wl, model_flux0, model_flux3, abs_err=True)
In [95]:
plot_diff_spec(wl, model_flux0, model_flux2, model_flux3, abs_err=True)
Measured in order 23
It seems like at this stage, with the Mg b lines, that metallicity contributes to the greatest change in the spectrum.
So this is interesting, a 1.5% change in flux level corresponds to a 0.05 dex change in metallicity.
Measured in order 23
Measured in order 23
Measured in order 41
The previous difference spectra were determined using our most accurate implementation of the spectrum modifier yet, that means using the full-resolution PHOENIX grids and k=5 polynomials, such that the spectra are accurate to within floating point precision. Since we have interpolated near a grid point, this also means that any linear approximation to the spline interpolation should be accurate enough.
To test any new approximate scheme for generating spectra, we want to see whether the residual spectrum will be bounded by the smallest spectrum in the {temp, logg, Z} trio, meaning that any error in the approximation will at most contribute to 10% error in any of {temp, logg, Z}. We can increase this threshold if it makes sense.
In [ ]: